Scattering of surface waves by a vertical truncated structured cylinder

This paper describes the solution to the problem of scattering of plane incident waves on water of constant depth by a bottom mounted circular cylinder, extending partially through the depth, which has an internal structure comprised of closely spaced thin vertical barriers between which fluid is allowed to flow. The problem is solved under full depth-dependent linearized water wave theory using an effective medium equation to describe the fluid motion in cylinder and effective boundary conditions to match that flow to the fluid region outside the cylinder. The interest in this problem lies in the development of novel solution methods for fully three-dimensional water wave interaction with bathymetric plate arrays. Results computed using this theory are compared with a shallow water approximation based on the recent work of Marangos & Porter (2021 Shallow water theory for structured bathymetry. Proc. R. Soc. A 477, 20210421.) and with accurate computations of an exact representation of the geometry using a discrete set of plates. Other results highlight the resonant directional lensing effects of this type of cylindrical plate array device.

This paper describes the solution to the problem of scattering of plane incident waves on water of constant depth by a bottom mounted circular cylinder, extending partially through the depth, which has an internal structure comprised of closely spaced thin vertical barriers between which fluid is allowed to flow. The problem is solved under full depth-dependent linearized water wave theory using an effective medium equation to describe the fluid motion in cylinder and effective boundary conditions to match that flow to the fluid region outside the cylinder. The interest in this problem lies in the development of novel solution methods for fully three-dimensional water wave interaction with bathymetric plate arrays. Results computed using this theory are compared with a shallow water approximation based on the recent work of Marangos & Porter (2021 Shallow water theory for structured bathymetry. Proc. R. Soc. A 477, 20210421.) and with accurate computations of an exact representation of the geometry using a discrete set of plates. Other results highlight the resonant directional lensing effects of this type of cylindrical plate array device.

Introduction
Closely spaced arrays of thin plates extending vertically from the bottom of an ideal fluid with a free surface have been used recently by a number of authors to produce unusual effects on water waves that are inaccessible using conventional changes in depth, notably the ability to negatively refract waves. Contributions include Berraquero et al. [1], Maurel et al. [2,3] and Marangos & Porter [4], all of whom assumed long wavelengths compared to the depth and subsequently developed depth-averaged models to describe the wave propagation over this plate-array structured bathymetry. Each of those models produces the same governing twodimensional wave equation in which the conventional fluid depth is replaced by a diagonal rank two tensor but whose elements differ depending on further assumptions made about the spacing of gaps within the plate array. The two diagonal tensor entries determine the phase velocities of surface waves propagating in different directions over the bathymetry and this anisotropy leads directly to negative refractive effects. It is also an important ingredient in bathymetric cloaking of vertical cylinders [5].
The model of Marangos & Porter [4] applies to close spacing between the plates where tensor elements are expressed explicitly in terms of the distances from the surface to the top of the array of plates and to the base of the fluid. This model arises from homogenization which exploits the contrast in the lengthscales in the problem and is adopted from Porter [6] and Zheng et al. [7]. In these two problems, the plate arrays extend throughout the fluid depth allowing the exact depth dependence of the fluid to be factorized from the solution without depth averaging and is therefore not restricted to long wavelengths. This factorization implies the water wave equations are analogous to the equations that govern two-dimensional acoustics and polarized electromagnetics.
Whilst Porter's [6] work concentrated on refraction across planar interfaces, Zheng et al. [7] considered the effect on waves of plate arrays confined within a circular cylinder. The small gaps between the plates provide an environment for resonance and slow wave propagation across the cylinder giving it the ability to concentrate wave energy which can either be harnessed/dissipated via a local damping mechanism or redirected to create a water wave lensing device. Moreover, when the incident wave direction is aligned with the plates, the cylinder becomes completely transparent to waves, making the cylinder an interesting prospect for marine energy harvesting since it can be rotated to protect itself when necessary whilst possessing the potential to capture significant amounts of energy when engaged against incident waves.
In this paper, we consider the effect that truncating the vertical extent of cylinder considered by Zheng et al. [7] has on the scattering of incident waves. Specifically, we have assumed the cylinder extends upwards from the base of the fluid to a constant level below the surface. Although this particular configuration is not a candidate for wave energy capture, it does allow results to be compared to the shallow water model of Marangos & Porter [4]. This serves an important purpose and is one of the motivations for the present work: calibrating the conditions under which the much simpler and more widely used shallow water approximation can be used to accurately determine scattering by bathymetric plate array devices.
The structural uniformity in the depth exploited by Zheng et al. [7] is absent here and this introduces additional mathematical challenges not previously encountered. It is the description of this novel and bespoke solution process for the full depth-dependent model on which the emphasis of this work is placed. We regard this as an important prototype problem for developing the solution methods for more complex problems involving wave energy harvesting by multiscale devices which involve mechanical damping components. A particular motivation for considering the bottom mounted submerged truncated cylinder as opposed to cylinders intersecting the surface is to avoid problems associated with undamped resonance for sufficiently high frequencies as encountered by Zheng et al. [7].
Separation solutions are employed but, within the cylindrical region defined by the cylinder radius, the field equation satisfied by the velocity potential above and below the submerged level of the top of the cylinder switches from the three-dimensional Laplace equation to a reduced

Description of the problem
We work with a mixture of three-dimensional Cartesian and cylindrical coordinates with z = 0 coinciding with the mean free surface of the fluid, which rests above a horizontal bed at z = −h.
Otherwise (x, y) = (r cos θ , r sin θ ) lie in the horizontal plane and a structured cylinder is enclosed The internal structure of the cylinder comprised thin vertical barriers which extend to the boundaries of the cylinder and are separated from one another by small uniform gaps through which the fluid is allowed to flow (figure 1). Without loss of generality, the plates are aligned with the x-axis since we allow a plane wave to be incident from infinity at an arbitrary angle, β. The incident wave is of angular frequency ω and wavenumber k and described, under classical water wave theory, by the velocity potential φ(r, θ, z) = e ikr cos(θ−β) ψ 0 (z), (2.1) (a time-harmonic dependence e −iωτ , in which τ is time, is suppressed hereafter) where In the above, k is related to ω and the depth, h, by the usual dispersion relation In (2.2) N 0 is a normalization factor defined by Thus φ inc is a solution of the governing equation in the fluid, (∇ 2 being the three-dimensional Laplacian) satisfying the bottom boundary condition, 6) and the combined linearized dynamic and kinematic boundary condition,  wavelength and the dimensions of the cylinder allow us to replace the microstructure by the effective field equation (formally derived in Porter [6]) where we have writtenφ(x, y, z) ≡ φ(r, θ , z). The zero normal flow conditions on the individual elements of the array are taken account of in the derivation of this reduced Laplace's equation. It is easy to interpret (2.8) as acting to restrict the fluid motion within the gaps in the cylinder to the x − z plane only. We must also consider how the field within the effective medium governing by (2.8) connects to the fluid outside the cylinder. This can be done formally, but amounts (at leading order in the small parameter on which the derivation of (2.8) is based) to matching the local pressures and local fluxes across the boundary of the cylinder. Thus, on z = −d, r < a, 0 < θ ≤ 2π we havẽ φ(r cos θ , r sin θ, −d) = φ(r, θ , −d) andφ z (r cos θ, r sin θ, −d) = φ z (r, θ, −d). (2.9) Over the curved surface of the cylinder, r = a, −h < z < −d, 0 < θ ≤ 2π , the matching conditions areφ (a cos θ , a sin θ, z) = φ(a, θ , z) and cosθφ x (a cos θ, a sin θ, z) = φ r (a, θ, z). (2.10) The term cos θ is geometric and arises from the conservation of mass flux across local triangular matching regions between the channel aligned with the x-axis and the radial flow into r > a.
(e.g. [7]). The only other constraint on φ is that it satisfies a standard radiation condition to ensure everything apart from the incident wave is radiating energy towards infinity.

Solution
The solution of the problem will be expressed using separation of variables inside and outside the cylindrical surface r = a and subsequently completed by matching across r = a appropriately. In r > a, we follow the standard method of expanding the incident wave into polar coordinates using the Jacobi-Anger expansion, in terms of Bessel functions, J n , and then the total potential, incorporating waves outgoing to infinity, is written in its most general form as i n e inθ a n,0 H n (kr)ψ 0 (z) + ∞ m=1 a n,m K n (k m r)ψ m (z) , (3.2) where a n,m are Fourier-Bessel expansion coefficients, to be determined, H n ≡ H (1) n are first-kind Hankel functions and K n are second-kind modified Bessel functions. We have defined depth eigenfunctions by where k = ik m are roots of (2.3) for m = 1, 2, . . ., lying on the positive imaginary k-axis: that is, defined by the positive real roots of K = −k m tan k m h such that (m − 1 2 )π < k m h < mπ for m = 1, 2, . . .. In doing so the orthogonality condition holds for n, m = 0, 1, 2, . . ., where δ n,m denotes the Kronecker delta, which is 1 if n = m, and 0 otherwise. The notation has been extended to include (2.2) via k 0 ≡ −ik. By contrast, determining the solution in r < a is more complicated on account of there being two distinct domains above and within the structured cylinder within which the governing equations differ, their solutions connected by the conditions (2.9).
We write the solution in r < a satisfying (2.5) in −d < z < 0 and (2.8) in −h < z < −d in its most general form, being a superposition over all possible wavenumbers and wave angles, thus φ(r, θ , z) = ∞ q=0 π −π B q (t)e iμ q (t)r cos(θ−t) Z q (z, t) dt, (3.5) where B q (t) are undetermined functions. The depth variation is defined in a piecewise fashion designed to satisfy (2.6), (2.7) and the matching conditions (2.9) on z = −d, by such that μ = μ q (t) are solutions of For every value of t in [−π , π ) there exist an infinite number of discrete roots of (3.8) labelled μ = μ q (t) where q = 0, 1, 2, . . .. It is shown in appendix A that these consist of two real roots labelled μ = ±μ 0 (t) and an infinite sequence of roots lying on the imaginary axis at μ = ±μ q (t), q = 1, 2, . . .. It is also shown in appendix A that there are no roots lying off the real and imaginary axes. We need only include the single positive real root and the sequence of roots that lie on the positive imaginary axis in (3.5) since these contribute to propagating and evanescent waves heading in the direction t, a variable which is integrated over all angles, −π ≤ t < π. Using the governing equations satisfied by Z q in −d < z < 0 and −h < z < −d, it can be readily shown that these functions satisfy the generalized orthogonality condition on account of Z q being real; C q (t) can be determined explicitly. The identity (3.9) is not used beyond this point, but could be useful in other problems. It is helpful to employ the Jacobi-Anger expansion of the plane wave function so that we have, in place of (3.5), the general expansion Now φ is continuous across r = a for all −h < z < 0, 0 < θ ≤ 2π . So from (3.2) and (3.11), matching Fourier modes in θ , multiplying by each of the depth eigenfunctions ψ m (z) and integrating over the depth gives and a n, and where which can be determined explicitly (see appendix B). The matching of fluxes across r = a is more complicated since the condition changes across 17) and this allows us to write These two expressions are matched over their respective intervals of depth to φ r | r=a + calculated from (3.2) which results in and k m a n, We can now eliminate a n,0 between (3.12) and (3.19), the resulting equation holding for −∞ < n < ∞, and a n,m from between (3.13) and (3.20), resulting in equations for −∞ < n < ∞ and for m = 1, 2, . . .. Combining these results gives rise to the system of equations (after using Abramowitz & Stegun ( [11], §9.1.6)) and

(a) Numerical approximation
We have to solve (3.21) for the functions B q (t), −π ≤ t < π for q = 0, 1, . . .. Since M q,n,m (t + 2π ) = M q,n,m (t) is a continuous smooth function we assume we can write The exponential factor involving the argument μ q (t)a is introduced to suppress the exponential behaviour of the functions J n (μ q (t)a) for q ≥ 1 when μ q (t) is imaginary, essential in avoiding numerical solutions becoming dominated by rounding errors. In particular, the NAG libraries, In the above, the assumed expansion of M q,n,m (t) in (3.25) has implied the expansion of B q (t), although an alternative approach would have been to expand B q (t) in terms of any basis whose elements are periodic in t with period 2π and determine a discrete system of equations that result. One final computational issue is that the equations nested within (3.18) can be multiplied by K n (k m a) for m = 0 (and H n (ka) when m = 0) and this suppresses a second potential source of exponential behaviour of O(e k m a ) from the elements M p,q,n,m .

(b) The far field diffraction coefficient
As kr → ∞ we have from (3.2) and introducing the large argument asymptotics of the Hankel function, that The scattering cross section, representing the total energy in circular waves diffracted by the cylinder, is defined as where Re denotes the real part of a complex number, and the last equality follows by the 'optical theorem' (Maruo [12] or see Mei [13], eqn (6.33)). Using (3.31) in (3.32) we have This latter relation is particularly useful for assessing the accuracy of solutions as it provides two independent calculations, σ 1 and σ 2 , say, of the same quantity.

(c) Forces
The net force acting on the structured cylinder is sum over all of the vertical barriers of the differential pressure acting over each barrier. It is straightforward to determine the effective medium limit of this discrete description which results in the expression after using a standard recurrence relation for derivatives of Bessel functions expressed in the form (xJ 1 (x)) = xJ 0 (x). Upon using the expansion (3.29) we have For the purposes of presentation, we plot the dimensionless quantity F y = F y /F cyl where F cyl is the force in the y-direction on a solid circular cylinder of radius a extending through the depth (e.g. [14]) subject to waves incident at an angle β to the positive x-axis and given by

Shallow water theory
The truncated structured cylinder protrudes from the bed and acts as a bathymetric metamaterial which can be considered on the assumption that incident waves are long compared to the fluid depth using the shallow water approximation of Marangos & Porter [4]. The approximation is therefore designed to work under the assumptions kh 1 and a/h 1 and, consequently, the field variable φ 0 (r, θ ) ≈ φ(r, θ , 0), proportional to the surface elevation, is written in r > a where k 2 h = K defines the wavenumber as the long wave limit of the dispersion relation (2.3). Inside r < a the governing shallow water equation for the metamaterial depth is, according to Marangos & Porter [4], is the two-dimensional gradient and the ∇ xy · is the Cartesian divergence operator. In (4.2) is a Cartesian tensor. The full depth of the fluid in r < a is denoted by D and this need not now be the same as h, the depth in r > a.
While the Cartesian description of the field equation is sufficient to generate a general solution within r < a, the boundary conditions on r = a require us to make a transformation into polar coordinates.
where ∇ rθ ≡ (∂ r , r −1 ∂ θ ) is the gradient and ∇ rθ · ≡ (r −1 ∂ r r, r −1 ∂ θ ) the divergence in polars while is the transformation of the Cartesian depth tensor into polars. Thus we have The matching conditions at r = a are: (i) that φ is continuous; and (ii) that the depth-averaged flux is continuous. The latter condition is equivalent to (since the flux vector with components in the radial and angular directions respectively is h rθ ∇ rθ φ 0 in transformed coordinates).

(b) Solution
This is the analogue of the method used for the full depth-dependent model without the complication of the depth variation. The general solution of (4.2) in r < a can be written as a superposition over plane waves travelling in all directions where, in order that the governing equation (4.2) be satisfied, It can be confirmed that (4.8) satisfies (4.6) also. Expanding the complex exponential as a series over Bessel functions (3.10) gives us We can apply the matching conditions relatively easily now. The matching of φ 0 (a, θ) for 0 ≤ θ < 2π results in J n (ka)e −inβ + a n H n (ka) = The flux condition derived in polars in (4.7) must be used to generate a second relation between the coefficients a n and b n and is clearly more complicated than the first condition. We find, after some algebra that kh(J n (ka)e −inβ + a n H n (ka)) = π −π Bessel function recurrence relations (J n−1 (x) − J n+1 (x) = 2J n (x), and J n−1 (x) + J n+1 (x) = 2nJ n (x)/x) provide simplification to kh(J n (ka)e −inβ + a n H n (ka)) We eliminate a n from between (4.12) and (4.13) to get π −π B(t)M n (t) dt = G n,0 , (4.14) where G n,0 is defined by (3.22) and where . This approach follows the one we used for the fully depth-dependent formulation in which we have not sought to approximate B(t) directly. The choice of expanding M n (t) in a Fourier basis implies that the solution of (4.18) encodes the Fourier coefficients of B(t) since it follows from (4.19) that (4.20) The two expressions presented in (3.33) can be used as approximations σ 1 and σ 2 to the scattering cross section, σ , with a n,0 replaced by a n .
The horizontal force on the truncated structured cylinder under the shallow water assumptions is This should be normalized by the shallow water limit of (3.39) in which the value of the square brackets is set to unity.

Results
In both the full depth-dependent treatment of the problem and the much simpler shallow water approximation, there are a number of numerical parameters whose influence on the accuracy of computations needs to be considered. Thus, range of values of n and p representing the angular variation are reduced to −N ≤ n, p ≤ N and the range of values of m and q, representing the vertical variation (or the number of evanescent modes retained), are reduced to 0 ≤ m, q ≤ M. Apart from these parameters, integrals require numerical approximation. The integrals defined over −π ≤ t < π are first arranged as integrals over 0 ≤ t < π/2 and then treated using a non-adaptive Gaussian quadrature which is refined to ensure that results are free from quadrature errors up to the sixth decimal place. It has been confirmed that in the special cases of β = 0, π there is no scattering for any a/h, d/h (in addition to D/h = 1 in the case of shallow water) and this trivial result is insensitive to truncation parameters N, M. There is no value in presenting numerical results relating to such cases. Table 1 presents tabulated results that are designed to illustrate the convergence of the numerical scheme with N and M in a non-trivial case. The results are typical of those found for other values of parameters a/h, d/h and β provided either d/h is not too small or a/h is not excessively large. The task is made difficult by the sensitivity in computational stability to increasing values of N and M and neither can become too large without the accuracy of the computations becoming compromised by numerical rounding errors. In particular, increasing N beyond the value of 12 can easily result in numerical errors. Identifying the source of these numerical errors has been difficult. The numerical solution requires the computation of Bessel functions (NAG libraries) of both large argument and large order, the numerical integration of functions which oscillate with increasing frequency with N and have increasingly abrupt changes in μ q (t) as q increases and the inversion of the complex matrix of increasingly large dimension (2N + 1) × (M + 1). On the other hand, we can see from table 1 that numerical results converge sufficiently well for presented curves and surfaces to be graphically accurate (i.e. two or three decimal place accuracy) with relatively small values of N and M. The graphical results produced in the paper have used values of N and M within the range of values presented in table 1. Generally, for larger values of ka we require larger N to represent higher frequency diffracted wave effects but generally smaller M. Conversely, when ka is small, N can take small values but M should be larger as the structure influences the fluid motion at larger depths when subject to longer waves. Figure 2 shows the variation of σ , the scattering cross section, and | F y |, the magnitude of the dimensionless force, computed using full depth-dependent theory and shallow water theory. In the upper two subplots a/h = 4, d/h = 0.5 and the height of the cylinder one-sixteenth its diameter and in the lower two subplots a/h = 0.5, d/h = 0.1 and the height is 90% of the cylinder diameter. We see that the shallow water approximation is in good agreement with the full linear theory for smaller values of kh as expected. The shallow water approximation is not designed to work with large and abrupt changes in depth and will therefore tend to work better for short wide cylinders a/h 1 rather than tall narrow cylinders. The oscillatory behaviour of the force as a function of kh in figure 2b is attributed to multiple interference effects of the waves over the top of the cylinder. We note that the force on the structure is of the same order of magnitude as that on an equivalent solid cylinder of the same size.
Further comparisons are made between full depth-dependent theory and shallow water theory in figure 3 for kh = 1, beyond the wavenumber regime, kh 1, where shallow water approximation is designed to work. Nevertheless, there is a good qualitative agreement although there are observable differences in the wave elevation particularly in the large amplitudes across the top of the cylinder in the case of d/h = 0.1. These large amplitudes are not simply a shoaling effect due to the reduced depth of water; in both plots, we see the signature of the resonance which exists when plates extend fully throughout the depth, as reported by Zheng et al. [7].     This near-resonant behaviour is amplified as d/h decreases and is therefore more prominent in figure 3a.
The accuracy of the effective medium model developed in this paper is tested in figure 4 in a comparison with a computation using the boundary element method of Liang et al. [15] for an arrangement of discrete thin vertical plates. Boundary element method computations were performed at different wavenumbers and for arrangements of 10, 20 and 40 plates and showed convergence to the effective medium model with increasing numbers of plates. Just one example has been used here for illustration in figure 4 for 20 plates and for a cylinder extending through 80% of the depth. It can be seen that there is visibly almost perfect agreement between the surface plots both inside and away from the cylinder, with the effective medium results predicting slightly more amplification than the discrete computation.
The lensing effect of the plate array cylinder is highlighted by the plots in figure 5 showing the magnitude of the diffraction coefficient |A(θ )|, which measures the wave amplitude of circular waves propagating outwards in the direction θ against the incident wave direction β. In figure 5   the lowest value of kh = 1 shown there is an overall reduction in wave scattering but with increased back-scattering into θ ∈ (−180 • , 0 • ), the diffracted becoming increasingly symmetric about θ = 0. The diffraction in the long wavelength limit tends to a dipole and this is confirmed by computations which show that the size of the coefficients, a n,0 , contributing to outgoing wave propagation, are increasingly dominated by a ±1,0 as kh → 0.
Returning to figure 5 for kh ≈ 1 − 2 see diffracted amplitudes distributed mainly about a maximum close to θ = 90 • across a wide range of values of β with decreasing amounts of backscatter into θ ∈ (−180 • , 0 • ) as kh increases. As kh increases there is a tendency for incident waves  to be scattered forwards in the same direction as the incident wave. In all cases, as β approaches 90 • we find the forward scattering at θ = 90 • reaches a maximum and the magnitude of this scattering increases with kh. The features described above are largely in common with the plate array cylinder extending through the depth described by Zheng et al. [7]. Thus, we have shown that the submerged truncated plate array cylinder can be used as a local resonator with the capacity to accumulate wave energy which is transmitted away from the cylinder in a direction perpendicular to the plate array. Finally, in figure 6 we showcase the ability of the the shallow water approximation to deal with the case of a cylindrical hollow which extends in r < a to a level D = 2h below the depth h in r > a. The hollow is filled with a plate array up to the level d = h. When β = 90 • waves of any frequency are transparent to this plate array 'pit'. In figure 6 where β = 0 • there is scattering creating a quiet zone in the lee of the 'pit', which is more a function of the lower depth in the pit than due to the presence of the plate array within the pit.

Conclusion
In this paper, we have considered wave scattering by a porous vertical bottom-mounted cylinder extending part of the way through the depth. The internal structure of the cylinder comprised closely spaced thin vertical barriers whose effect is to confine the motion of the fluid flow in the narrow gaps between the barriers. This allows us to consider that the field inside the truncated cylinder is governed by an effective medium equation with effective matching conditions holding on its boundary. The main focus of the paper has been to develop a solution to the full depthdependent potential flow problem which exploits this effective medium description of the cylinder. Results have been tested against the shallow water approximation of Marangos & Porter [4] and good agreement is reached for sufficiently long wavelengths as expected. Results are also compared to a boundary element computation of an exact geometrical description involving a finite number of discrete barriers based on the work of Liang et al. [15], again showing excellent agreement. Near-resonant wave trapping is promoted within and above the truncated cylinders for sufficiently large wavenumbers. The truncation is important since it implies our solutions do not suffer from undamped resonances reported by Zheng et al. [7] for plate array cylinders extending fully throughout the depth. Like the work of Zheng et al. [7], we have shown that a by-product of this near-resonant trapping is directional scattering, or lensing, of wave energy perpendicular to the alignment of the plate arrays. New results show that there is no scattering of wave energy in directions parallel to the plate array and that, for long wavelengths (compared to the cylinder radius), the dominant scattering pattern is dipolar. This latter result, which contrasts with the monopolar scattering for small solid cylinders, highlights that small plate array cylinders could be used as a physically realizable element for dipolar scattering in, for example, metamaterial design.
The topic of the current paper and the development of the mathematical methods herein has been considered with a particular, more practical, extension in mind. Work now being considered by the authors involves developing a cylindrical wave energy converter involving a truncated plate array cylinder which intersects the free surface.
Data accessibility. Numerical code and data files for the results in this article are available to download from https://people.maths.bris.ac.uk/~marp/abstracts/metatrunccyl.html. This appendix concerns the roots of (3.8) for −π ≤ t < π. First, we note that the roots are symmetric about t = 0 and t = π/2 and therefore need only be found for 0 ≤ t ≤ π/2. For t = π/2 the roots of (3.8) coincide with those of K = μ tanh μd, while for t = 0 the roots coincide with those of K = μ tanh μh, the dispersion relation for waves propagating in fluid of depth d and h, respectively. Thus, in these two cases, the locations of roots are well understood (e.g. [16]): there are two real roots and an infinite sequence of imaginary roots, symmetric about real and imaginary axes in the complex plane; there are no complex roots.
The remainder of this appendix concerns the roots of (3.8) for values of t ∈ (0, π/2) away from these special cases. This is important in developing numerical schemes to ensure all possible roots are captured and can be computed efficiently. We approach this with care since, for example, the dispersion relation corresponding to waves in the presence of a thin floating elastic plate on water of finite depth gives rise to roots which, under certain conditions, move away from the imaginary axes into the complex plane as physical parameters vary. In the elastic plate case, this behaviour is associated with the coalescing of pairs of roots on the imaginary axis [17]. However, in the analysis below we will show that the roots of (3.8) are located on either the real or imaginary axes for all t ∈ (0, π/2).
This result provides us with our first (and simplest) numerical method for determining roots. Since we are required to compute integrals over 0 ≤ t < π/2 we determine the roots of the water wave dispersion relation for t = 0 and then track the roots as t is varied across the integration range, ensuring at t = π/2 the expected roots are accounted for and none have been 'lost'.

(b) Imaginary roots
To consider roots of (3.8) lying on the imaginary μ axis let μ = iμ whereμ is real such that (3.8) which can be writtenL (μ) =R(μ), We note thatL andR are odd functions and so ifμ > 0 is a root then so is −μ.
These features are illustrated in figure 8 whence it can be concluded that there are roots of (3.8) between each asymptote ofL andR countingμ = 0 as the first asymptote when Kd ≥ 1, or the first  asymptote of eitherR orL when Kd < 1. Thus discrete roots exist within intervals that are known explicitly and this allows us to solve for roots using, for example, bisection. This provides us with an alternative numerical scheme to the one outlined previously. We have checked that both methods give the same results. An example of how the roots vary as a function of t is provided in figure 9.
Then F 1 and F 2 are meromorphic functions. We will use Rouché's Theorem which states that the number of zeros of F = F 1 + F 2 inside a closed contour, C, in the complex plane is equal to the number of zeros of F 1 inside C provided |F 1 | > |F 2 | for all points on C.
The choice of the specific distance at which C m (ρ) cuts the imaginary axis is key to the result. Thus it can be confirmed, again most easily by graphical considerations, that between −imπ/[(h − d) cos t] and +imπ/[(h − d) cos t] on the imaginary axis there are as many zeros of F = F 1 + F 2 as there are zeros of F 1 = R 1 L 2 . For example, referring to figure 8, we choose m = 2 in which C m (ρ) cuts the imaginary axis at the second positive zero ofR(μ) orμd ≈ 8.9. We can see there are four roots of F, denoted by circles, lying on the positive imaginary axis in the interval 0 ≤μd 8.9. In the same interval, there are two zeros ofR, coinciding with the zeros of R 1 , and two asymptotes ofL, which coincide with the zeros of L 2 , and thus four zeros of F 1 . Also, there are exactly two roots of F symmetric placed on the real axis, the same as the number of roots of F 1 on the real axis (L 2 has no zeros but R 1 has two zeros symmetrically placed about the origin; see figure 7).
It can also be shown easily that |L 2 | > |L 1 | on C m (ρ) and (less easily analytically, but also confirmed numerically) that |R 1 | > |R 2 | on C m (ρ) as ρ → ∞ for each m. Thus, |F 1 | > |F 2 | on C m (ρ) as ρ → ∞ as required by Rouché's Theorem. Finally, since we know that the zeros of R 1 only lie on the imaginary axis and the zeros of L 2 only lie on the real and imaginary axes (being roots of the water wave dispersion equation for water depth d), it must follow that the zeros of F also only lie on the real and imaginary axes.